import numpy as np
import matplotlib.pyplot as plt

plt.rcParams['font.sans-serif'] = ['SimHei']  # 用黑体显示中文
plt.rcParams['axes.unicode_minus'] = False    # 正常显示负号

# 实验一：解析函数和调和函数的验证 (简化版)
print("实验一：解析函数和调和函数的验证")
print("="*40)

# 定义复函数 f(z) = z^2
def f(z):
    return z**2

# 验证复函数的解析性（通过Cauchy-Riemann方程）
def check_analyticity(func, x, y, h=1e-6):
    z = x + 1j*y
    # u_x = (u(x+h, y) - u(x, y))/h, v_y ≈ (v(x, y+h) - v(x, y))/h
    u_x = (np.real(func((x+h) + 1j*y)) - np.real(func(z))) / h
    v_y = (np.imag(func(x + 1j*(y+h))) - np.imag(func(z))) / h
    # u_y = (u(x, y+h) - u(x, y))/h, v_x ≈ (v(x+h, y) - v(x, y))/h
    u_y = (np.real(func(x + 1j*(y+h))) - np.real(func(z))) / h
    v_x = (np.imag(func((x+h) + 1j*y)) - np.imag(func(z))) / h
    # Cauchy-Riemann方程: u_x = v_y, u_y = -v_x
    return abs(u_x - v_y) < 1e-4 and abs(u_y + v_x) < 1e-4

# 验证调和函数性质
def check_harmonic(u, x, y, h=1e-4):
    # u_xx = (u(x+h, y) - 2*u(x, y) + u(x-h, y))/h^2
    u_xx = (u(x+h, y) - 2*u(x, y) + u(x-h, y)) / (h**2)
    # u_yy = (u(x, y+h) - 2*u(x, y) + u(x, y-h))/h^2
    u_yy = (u(x, y+h) - 2*u(x, y) + u(x, y-h)) / (h**2)
    # 拉普拉斯方程: Δu = u_xx + u_yy = 0
    return abs(u_xx + u_yy) < 1e-3

# 对于 f(z) = z^2 = (x+iy)^2 = x^2 - y^2 + 2ixy
def u(x, y): return x**2 - y**2  # 实部
def v(x, y): return 2*x*y         # 虚部

# 验证解析性
print("1. 验证 $f(z) = z^2$ 的解析性:")
point_x, point_y = 1.0, 1.0
is_analytic = check_analyticity(f, point_x, point_y)
print(f"  函数在点 ({point_x}, {point_y}) 解析: {is_analytic}")

# 验证调和性
print("\n2. 验证调和函数:")
is_harmonic_u = check_harmonic(u, point_x, point_y)
is_harmonic_v = check_harmonic(v, point_x, point_y)
print(f"  $u(x,y) = x^2-y^2$ 是调和函数: {is_harmonic_u}")
print(f"  $v(x,y) = 2xy$ 是调和函数: {is_harmonic_v}")

# 可视化调和函数
print("\n3. 可视化调和函数 $u(x, y) = x^2 - y^2$:")
x = np.linspace(-2, 2, 50)
y = np.linspace(-2, 2, 50)
X, Y = np.meshgrid(x, y)
Z = u(X, Y)

fig, ax = plt.subplots(1, 1, figsize=(6, 5))
contour = ax.contour(X, Y, Z, levels=15)
ax.clabel(contour, inline=True, fontsize=7)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_title('调和函数 $u(x, y) = x^2 - y^2$')
ax.grid(True)
plt.tight_layout()
plt.show()